Quantum paramagnets on the honeycomb lattice and field-induced Neel order: 

Possible application to Bi3Mn40i2(N03) 
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Motivated by recent experiments on the spin-3/2 frustrated bilayer honeycomb antiferromagnet 
Bi3Mn40i2(N03), we study the spin-S Heisenberg model on the honeycomb lattice with various 
additional exchange interactions which frustrate Neel order. Using spin wave theory, exact diago- 
nalization, and bond operator theory, we consider the effects of (i) second-neighbor exchange, (ii) 
biquadratic exchange for S = 3/2 which leads to an AKLT valence bond solid, and (iii) bilayer 
coupling which leads to an interlayer dimer solid. We show that the resulting paramagnetic states 
undergo a transition to Neel order beyond a critical magnetic field. We discuss experimental impli- 
cations for Bi3Mn40i2(N03). 



The interplay of quantum mechanics and frustrated 
interactions in quantum magnets leads to a variety of 
remarkable phases including spin liquid Mott insulators, 
valence bond crystals, and Bose-Einstein condensates of 
magnons [ij. Experiments on Bi3Mn40i2(N03) indi- 
cate that this material is a possible new candidate for a 
quantum spin liquid [2,] . The octahedral crystal field, to- 
gether with strong Hund's coupling, leads to Heisenberg- 
like spin-3/2 moments on the Mn^"*" ions which form a 
bilayer honeycomb lattice. Despite the bipartite struc- 
ture, and a large antiferromagnetic Curie- Weiss constant 
Qcw ~ —257K, this system shows no magnetic order 
(or any other phase transition) down to T ~ IK Q . This 
observation hints at frustrating interactions which may 
lead to interesting paramagnetic ground states 

Recent neutron scattering experiments fll| on powder 
samples of Bi3Mn40i2(N03) in zero magnetic field indi- 
cate that there are short range spin correlations in this 
material, with some antiferromagnetic coupling between 
the two layers forming the bilayer, but negligible interac- 
tions between adjacent bilayers. Remarkably, a critical 
magnetic field, i?c ~ 6 Tesla, leads to sharp Bragg s pot s 
consistent with three dimensional (3D) Neel order 11 1. 
Motivated by the observation that the field required to in- 
duce Neel order appears to extrapolate to a nonzero value 
at T = 0, we propose that this system could exhibit a 
field tuned quantum phase transition into the Neel state. 
We flesh out this idea by studying various interactions 
which could frustrate the Neel order in this material. 

We first examine the possibility that the Neel order in 
Bi3Mn40i2(N03) is destroyed by a frustrating second- 
neighbor exchange ( J2) in addition to the dominant near- 
est neighbor term (Ji). Sidestepping the issue of what 
state results from such quantum melting, we study the 
magnetic field dependence of the critical J2 / Ji required 
to destroy the Neel order. Using spin- wave theory, we 
show that a nonzero magnetic field enhances the critical 
J2/J1, opening up a regime where applying a critical field 
to the non-Neel state yields long-range Neel order. 



We next explore other frustrating interactions which 
might kill Neel order and lead to novel quantum para- 
magnetic ground states on the honeycomb lattice. We 
focus here on two valence bond solid (VBS) states, 
which do not break any symmetries and are expected 
to show no thermal phase transitions as is the case with 
Bi3Mn40i2(N03). (i) Motivated by Bi3Mn40i2(N03), 
we study a generalized spin-3/2 model including bi- 
quadratic and bicubic spin interactions which permits an 
AfHeck-Kennedy-Lieb-Tasaki (AKLT) ground state [Tl - 
14| . Using exact diagonalization (ED) to compute the 
fidelity susceptibility |15|, we show that this model ex- 
hibits a direct Neel-AKLT transition. We also obtain 
the spin gap at the AKLT point, (ii) In view of the 
fact that Bi3Mn40i2(N03) consists of stacked bilayers, 
we use a spin-S generalization [l6| of the bond operator 
formalism [17] to show that a sufficiently strong bilayer 
coupling leads to an interlayer VBS state. Both valence 
bond solids, the interlayer VBS and the AKLT state, 
are shown to undergo a magnetic field induced quantum 
phase transition into a state which exhibits Neel order. 
Our results on the AKLT state are of broader interest 
given recent proposals to use this state as a universal 
quantum computation resource 18 1. 

Finally, we discuss possible experiments on 
Bi3Mn40i2(N03) which may help to discriminate 
between the various states we have studied, and to 
distinguish them from possible Z2 spin liquids ,9\. 

Second-neighbor exchange. — It has been suggested 
that the absence of Neel order in Bi3Mn40i2(N03) is 
linked to non- negligible further neighbor interactions [2} ■ 
We therefore study a minimal Hamiltonian, 



H 



S,-S^5f (1) 



where (.) and ((.)) denote nearest and next-nearest neigh- 
bor bonds respectively, and i? is a Zeeman field. Let us 
begin with a classical analysis valid for 5 = 00. When 
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J2=B=0, the ground state has colhnear Neel order. For 
J2 = and Bj^O, the spins in the Neel state start off 
in the plane perpendicular to the applied field and cant 
along the field direction until they are fully polarized for 
B > 6J1S. For B < 6J1S, the spin components trans- 
verse to the magnetic field have staggered Neel order for 
J2 < Ji/Q] this gives way to a one-parameter family of 
degenerate (canted) spirals for J2> Ji/60|- 

Incorporating quantum fluctuations is likely to lead to 
melting of Neel order even for J2 < Ji/6. Such fluctu- 
ations are also likely to completely suppress the classi- 
cal spiral order Using spin wave theory, we argue 
here that a small nonzero B enhances the stability of the 
Neel order compared to the zero field case, (i) For small 
nonzero B, spin canting leads to a small decrease, cx B^, 
in the classical staggered magnetization transverse to the 
field, (ii) On the other hand, one of the two magnon 
modes (labelled fij^) acquires a nonzero gap oc i? at the F- 
point as shown in Fig.[lja) (for S~3/2 with J2 — O.I5J1 
and B — O.5J1S). This suppresses quantum fluctuations. 
For B-^6JiS, the latter effect overwhelms the former, 
leading to enhanced stability of Neel order. 

To estimate the 'melting curve', we assume that the 
transverse spin components have Neel order along the 
S'jj-direction, and use a heuristic Lindemann-like crite- 
rion for melting: \/ {S^) — {Sx}"^ > a{Sx) where the ex- 
pectation values are evaluated in linear spin wave theory 
(see Supplementary Information for details). As shown 
in Fig[T](b) and its inset, quantum fluctuations at B = 
lead to melting of Neel order even for J2 < Ji/6 (i.e., 
before the classical destruction of Neel order). We set 
a=3 since this leads to a melting of Neel order for 5=1/2 
at J2 ~ O.O8J1, in agreement with a recent variational 
Monte Carlo study by Clark ei aZ 

For nonzero B, the 'melting point' moves towards 
larger J2 , leading to a window of J2 over which the quan- 
tum disordered liquid can undergo a field-induced phase 
transition to Neel order (see Fig. Hl^b)). This is consis- 
tent with recent neutron diffraction experiments [ll| on 
Bi3Mn40i2(N03). The window of J2 where such physics 
is operative appears to be small for S = 3/2: however, 
disorder effects would suppress the stiffness [l9j and may 
enhance this regime. We expect field induced Neel order 
even for 5 = 1/2 (see inset to Fig. [1Kb)). This can be 
verified by including a magnetic field in calculations re- 
ported in Ref.^. Our results also explain recent Monte 
Carlo simulations of the classical model with B^O [3] ; if 
J2=0.175Ji, as in the simulations, a nonzero B takes us 
closer to the melting curve, and may lead to the numeri- 
cally observed enhanced Neel correlations. Nevertheless, 
we expect that there will be no field-induced long-range 
Neel order for J2=0.175Ji in the classical model. 

Next-neighbor exchange thus provides a plau- 
sible explanation for the experimental data on 
Bi3Mn40i2(N03). We next turn to an exploration of 
other mechanisms which frustrate Neel order. 
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FIG. 1: (a) Dispersion of magnon modes Slj^ in the J1-J2 
model along depicted path in the Brillouin zone for J2=0.15Ji , 
S = 3/2 and B = 0.5JiS. (b) T = melting of Neel order for 
5=3/2 in the B-J2 plane (open triangles) obtained using a 
Lindemann-like criterion, \/ {S'}} — {Sx}^ ~3{Sx) ■ The region 
"???" is a quantum disordered state - possibly a valence bond 
solid or a quantum spin liquid. Arrow depicts path along 
which one obtains a field-induced transition to Neel order. 
Inset depicts a similar melting curve for 5*= 1/2. 

AKLT valence bond solid. — A particularly interest- 
ing spin-gapped ground state of a magnet with spin-S* 
atoms on a lattice of coordination number z—2S, is an 
AKLT valence bond state. Each spin-S* is viewed as 
being composed of 25" spin-1/2 moments symmetrized 
on-site, with each spi n-1/ 2 moment forming a singlet 
with one neighbor |12h14| . It was originally proposed 
as a realization of Haldane's prediction of a spin-gapped 
ground state in ID integer spin systems (20| . Assuming 
that the Mn''+ ions in Bi3Mn40i2(N03) mainly inter- 
act with the three neighboring spins in the same plane, 
this condition is satisfied with S = 3/2 and z = 3. 
The honeycomb lattice AKLT state has exponentially 
decaying spin correlations flS*!, and it is the exact, and 
unique, zero energy ground state of the parent Hamil- 



tonian Haklt = J2{^J)P^^■ 



Here Pj;^^^ denotes a pro- 
jector on to total spin-^ for a pair of spins on nearest 
neighbor sites (i,j). 
p.(3.) ^ 21. 



Sj, we find 
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We do not 
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11 I 243 rp , 116 Ti2 I 16 
2J 128 

have a microscopic basis, at this point, for such higher 
order exchange terms in Bi3Mn40i2(N03); but it is en- 
couraging to note that the coefficients of such terms are 
smaller than the leading Heisenberg interaction. 

We have investigated, using ED on system sizes N — 
12-18, the phase diagram of a generalized spin-3/2 model. 



Hq =(l-0)^S,-S, + (/Qi/AKLT, 



(2) 



which interpolates between a Heisenberg model (at Q=0) 
and gHAKLT (at Q=l). We set g=1440/243, so that the 
coefficient of Si-Sj is unity. For (5 = 0, our analysis of the 
finite size spectrum shows that the ground state energy 
Eg{N,S^°^), as a function of total spin 5*°', varies as 
5*°' (5'°* -1-1), in agreement with the expected Anderson 
tower for a Neel ordered state. It is consistent with ear- 
lier work showing Neel order even for spin-1/2 [s, 21, 2^ . 
To establish the Neel- AKLT transition as a function of 
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FIG. 2: (Color online) (a) Overlap P{Q\0) of the ground state 
at Q with the Neel state (Q=0) for various system sizes A'', 
showing its rapid drop around the Neel-AKLT transition, (b) 
Fidelity susceptibility xf{Q) versus Q for various system sizes 
A^, with the peak indicating the Neel-AKLT transition point 
Qc{N). Inset: Qc(A') versus 1/N, together with a fit Qc{N) = 
Q'^-\hN~^ (with a choice v ~ 0.7 assuming an 0(3) quantum 
phase transition in 2D) which leads to ~ 0.8. 

g, we Study overlaps P{Q\Q') = \{^ g{Q)Yl! g{Q'))\ of the 
ground state wave functions at Q and Q' . As shown in 
Fig Eta), the overlap P(Q|0), of the ground state wave- 
function at Q with the Neel state at Q' — Q, is nearly 
unity for (5<0.8, suggesting that the ground state in this 
regime has Neel character. For 0.8<(5<l-2, we observe 
a dramatic drop of P{Q\Q) for all system sizes, which 
indicates a Neel-AKLT quantum phase transition. To 
locate the transition more precisely, we compute the 'fi- 
delity susceptibility' xf{Q)= 2(1- P(Q|g (5))/(5^ with 
(5—^0, which measures the change of the wavefunction 
when Q [l5| . Vig^h) shows a plot of Xi^ (0) (with 

(5 = 0.005). We observe a peak in xf{Q) which indicates a 
phase transition; this peak shifts and grows sharper with 
increasing N . Assuming the thermodynamic transition 
is at Q^, and that the peak position QdN) satisfies the 
scaling relation (Qc(iV) - Q?") ~ iV"^/^'', with v w 0.7 
for an 0(3) quantum phase transition [il] corresponding 
to triplon condensation, we estimate the transition point 
g^wO.8. 

The spin gap A,{N) = Eg{N, S^°^ = l)~Eg{N, 5*°'=0) 
is plotted in FiglSl^a) for various Q as a function of 1/iV. 
Assuming a finite size scaling form As{N) = Af + b/N, 
we find a small value for for Q — 0.0,0.4, consistent 
with a gapless Neel state, while for Q — 0.9,1.0 there 
appears to be a robust spin gap as 1/A^ 0. At the 
AKLT point (Q=l), we estimate Af sa 0.6. 

Since the spin gap is finite for Q>Qc, we expect that 
applying a critical field Be oc A^ will lead to a phase tran- 
sition; the correlation functions of the 5'°* = 1 state at 
zero field will then reflect the correlations of the ground 
state for Bz>Bc. We plot, in Fig|3tb), the spin correla- 
tions on two maximally separated sites (for A^==16) as a 
function of Q, and make the following observations, (i) 
For 5**°* = 0, the ground state also has 5*°* = 0, and 
{Sx{i)Sx{j)) = {Sz{i)Sz{j)) due to spin rotational invari- 
ance. At long distance, the spin correlation is strong in 
the Neel phase, but drops rapidly to small values upon 
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FIG. 3: (Color online) (a) Spin gap, As(Af), versus l/A^ for 
various Q, with fits to the form As(A') = AJ° -I- 6/A'. The 
small values of for Q = 0.0, 0.4 are consistent with a gap- 
less Neel state. For Q = 0.9, 1.0, the data are consistent with 
a robust spin gap A^. (b) S^-spin correlations between dis- 
tant sites on the same {AA) and opposite (AB) sublattices 
for A''=16 system. The spin correlation is Neel-like (±) for 
all Q shown; in the spin gapped AKLT state at Q~l, it short 
ranged and weak in the S^°^ — ground state but it is strongly 
enhanced (see arrow) in the Sl°^ = 1 case. 

entering the AKLT state, (ii) In the 5**°* = 1 sector, 
(S'^(i)S'^(j)) ^ {Sxii)Sx{j)). Remarkably, in this sec- 
tor, as opposed to ^'"'^O, we find a strong enhancement 
of only transverse correlations {Sx{i)Sxij)) between dis- 
tant sites in the AKLT state; this finite-size result sug- 
gests that the AKLT state will undergo, beyond a critical 
field, a transition into a state with in-plane Neel order. 

Interlayer valence bond solid. — The Mn sites in a 
unit cell of Bi3Mn40i2(N03) form an AA stacked bilayer 
honeycomb lattice. If the interplane antiferromagnetic 
exchange Jc is strong compared to the in-plane nearest 
neighbor exchange Ji, adjacent spins on the two layers 
could dimerize and lead to loss of Neel order. To study 
this interlayer VBS, we begin from the limit Ji=0; this 
leads to the spectrum Ej = -Jc{S{S + 1) - j(j + 1)/2), 
with j = 0, 1, . . . , 25* denoting the total spin state of the 
dimer. Restricting attention to the low energy Hilbert 
space spanned by the singlet and the triplet states, we 
define generalized spin-S bond operators via: |s) = s^|0), 
and \a) = ^qIO); where |0) is the vacuum, and |a(= 
x,y,z)) are related to the ruj levels of the triplet by 
\z) = \mj = 0), \x) = {\mj = -l)-\mj = l))/\/2, and 
\y) — i{\m.j — — I) + \mj — 1)) / ^/2. Denoting the two spins 
constituting the dimer, by S^, with layer index £ = 0/1, 
we obtain [Ira 



S{S + 1] 



^ajSjt'pt^, (3) 



together with the constraint s^s + — 1 at each site. 

To treat the effect of Ji, we use bond operator mean 
field theory [l7| which yields a reasonably accurate phase 
diagram for the spin- 1/2 bilayer square lattice Heisen- 
berg model [i^. Assuming the singlets are condensed in 
the dimer solid, we replace s'^ = s = s, and incorporate 
a Lagrange multiplier in the Hamiltonian which enforces 
i^a^a) — 1 — on average. Let N be the number of spins 
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FIG. 4: Phase diagram of the S = 3/2 bilayer honeycomb 
model obtained using bond operator theory, and triplon dis- 
persion along depicted path in the Brillouin zone within the 
interlayer VBS state for Jc/Ji =3.8 (in units where Ji = l). 

in each honeycomb layer. We then obtain the Hamilto- 



Q,k>0 ^k. 



'2NC , describing the dy- 

Here *L = (4al4a2^-kal^-ka2) 

(with 1, 2 denoting the two sublattices in each layer) and 
the matrix Afk takes the form 



nian H = ^ 
namics of the triplets 



/A 

B 



Bk 
k ^k 
Bk 



V^k 





Bl 





Bk 



(4) 



neighboring site in the adjacent layer, the AKLT state 
would nucleate three 5"= 1/2 moments on neighboring 
sites in the same plane, while spinless impurities in spin 
gapped Z2 fractionalized spin liquids Q, do not generi- 
cally lead to local moments. Sharp dispersing triplet ex- 
citations expected in valence bond solids discussed here 
could be looked for using single-crystal inelastic neutron 
scattering; by contrast, a spin liquid may not possess such 
sharp modes. Specific heat experiments in a magnetic 
field could test for possible Bose-Einstein condensation 
of triplet excitations as a route to Neel order. 

Finally, dimer crystals with broken symmetry could 
also be candidate ground states in Bi3Mn40i2(N03) - 
in which case, disorder must be responsible for wiping 
out the thermal transition expected of such crystals. If 
BiMNO does support a valence bond solid ground state, 
disorder and Dzyaloshinskii-Moriya couplings (permitted 
by the bilayer structure) may be responsible for the ob- 
served nonzero low temperature susceptibility. This is an 
interesting direction for future research. 
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0906816 and DMR-0611562 (DNS). 



with Ak = Jc-fi-JcS{S+l) and Bk = l-fkJiS{S+l)s^ . 
Here we have defined 7k = 1 + c^''^ '' + e~*'^ (°+''), with 
unit vectors a = x,b = —x/2 + \/3y/2, and the constant 

c = -1is^-l)-l{Jc-^i-JcS{s+l))-^Jcs^sis+l). [1 

Diagonalizing this Hamiltonian leads to the ground state 
energy per spin Eg = 2 ^Ek>o(^k+ + Ck-) + C' where 
^u± = VMAu±2\Bu\). Setting dEg/ds^ = dEg/dfi = 0, 
we obtain the mean field values of s and which minimize 
the ground state energy subject to the constraint. Solv- 
ing these equations numerically, we find that the spin-S* 
interlayer VBS is a stable phase for Jc > Ji,[S] where 
J43/2] « 3.3J1 and J*[l/2] w O.66J1. Quantum Monte 
Carlo studies of this model would be valuable in firmly 
establishing the value of J* [5"] as a function of S. Fig. |4] 
shows the triplon dispersion of the S ^ 3/2 interlayer 
VBS state at Jc = 3.8Ji (in units where Ji = 1) along 
high symmetry cuts in the hexagonal Brillouin zone. For 
Jc < J*[3/2], or in the presence of a magnetic field which 
can close the spin gap in the VBS state for Jc > J* [3/2], 
the low energy triplon mode at the F-point condenses; its 
eigenvector is consistent with Neel order. 

Discussion. — Motivated by recent experiments on 
Bi3Mn40i2(N03), we have studied various honeycomb 
lattice spin models which support quantum paramagnetic 
ground states that undergo field-induced phase transi- 
tions to Neel order. Detailed NMR studies of isolated 
nonmagnetic impurities subsituted for Mn may help dis- 
tinguish between these states. The interlayer VBS would 
have an impurity induced 5 = 3/2 local moment on the 
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SPIN WAVE FLUCTUATIONS AROUND 
CANTED NEEL STATE 

The Hamiltonian, with a magnetic field and second 
neighbour exchange, is 

if = Ji ^ S, • S, + J2 51 • S, -BY^SI (SI) 

As discussed in the main body, when J2 < Ji/6 and 
B < 6J1S, Neel ordering is in the plane perpendicular 
to the magnetic field, but the spins also uniformly cant 
in the direction of applied field, to maximally gain Zee- 
man energy. The classical spin state can be characterized 
by Sr = S'(± cosx, 0, sinx) on the two sublattices. We 
now define new spin operators, denoted by Ti^a, via a 
sublattice-dependent local spin rotation 

sinx 


(-)"cosx 

where a = 1, 2, is a sublattice index and i sums over each 
unit cell. 

The ground state has all spins pointing towards the 
new local-^^ axis. To study spin wave fluctuations, we 
rewrite the T operators in terms of Holstein-Primakoff 
bosons as follows: 





The Hamiltonian can now be rewritten as iJ w Eci + 
Hqu ■ The classical energy Eci is proportional to S'^ , and 
the leading order quantum correction, Hqu, is of order S. 
We get the value of the canting angle x by demanding 
that terms of order S^/"^, which are linear in the boson 
operators, should vanish, which yields 



smx = 



B 



6Ji5 

The classical energy is given by 



Eci 



- Ji cos2x + - J2 



B 
'S 



smx- 



(S3) 



(S4) 



where N is the number of sites in the honeycomb lattice. 
We take the magnetic field B to be of order S, so that the 
Zeeman term —BSf is treated on the same level as the 
exchange terms Jij Sj-Sj . The leading quantum correction 



Bqu 
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k>0 



- Ji cos2x + 3J2 
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smx 



(S5) 



where 

with 
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; Hlr — 
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(S6) 



3Ji cos2x — 6J2, 

2 J2{C0S ka + COS fcb + COs(fca + fcfe)} - 

Ji7kSin2x= |i^k|e*''^ 
-Ji7kC0S^X> 



B 



■smx, 



where 7k is as defined in the main body. This Hamilto- 
nian can be diagonalized by a bosonic Bogoliubov trans- 
formation. The eigenvalues are given by 



f7± = v/(/k±|Fk|)2-|GkP. 



(S7) 



The transformation matrix is given by 

P : 



f/2x2 
t/2x2 



G2X2 '52x2 
'5'2x2 G2X2 



where 



t/2 
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71 



G2X2— 



cosh Q 
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'52x2— 



sinh( 







cosh(/) y'"-^'^^ \^ sinh( 
where the angles d and (f> are given by 

|Gk| 



(S8) 



,(S9) 



tanh 26* = 
tanh 2(/i = 



/k-|Fk|' 

-|Gk| 
/k + li^kl" 
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FIG. SI: Melting of Neel order for 5= 3/2 in the B-J2 plane 
for T = 0,0.2,0.4 (in units of Ji). To the left of the curve, 
there is stable canted Neel order. To the right, fluctuations 
melt the in-plane Neel order. 

The matrix P preserves the commutation relations of the 
bosonic operators and diagonalizes the Hamiltonian, giv- 
ing Pt HP = Diag{f^k 

The expectation value of spin moments can be calcu- 



lated in terms of this new basis. For example, the in- 
plane component of the spin is given by 



^E('5L=i)=(^ + V2)cos; 



cosx 



i 

E[cos26'{l+2nB(r2k,-)}+cos20{l-f2ni3(17k,+)}],(Sll) 



k>0 



where nsi-) denotes the Bose distribution function. For 
T 7^ 0, a small coupling along the third dimension is 
necessary to allow for a stable magnetically ordered state. 
We take this into account by imposing an infrared cutoff 
A which is of the order the interplane couling; modes 
with energy greater than A look like 2d spin waves. 

As J2 is increased, fluctuations around the Neel state 
increase. We expect the Neel state to melt when fluctua- 
tions become comparable to the magnitude of the ordered 
moment. As discussed in the main body, we choose our 
melting criterion to a/ (S^) — {Sx)'^=3{Sx) , which gives a 
critical Jij 3\ in agreement with recent variational Monte 
Carlo results for 5* = 1/2. The resulting melting curve for 
S = 3/2, at zero temperature and at small temperatures, 
is shown in Figlj 



